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Abstract 

Finding  polynomial  roots  rapidly  and  accurately  is  an  important  problem  in  many  areas  of 
signal  processing.  We  present  a  new  program  which  is  a  combination  of  Muller’s  and  Newton’s 
method.  We  use  the  former  for  computing  a  root  of  the  deflated  polynomial  which  is  a  good 
estimate  for  the  root  of  the  original  polynomial.  This  estimate  is  improved  by  applying  Newton’s 
method  to  the  original  polynomial.  Test  polynomials  up  to  the  degree  10000  show  the  superiority 
of  our  program  over  the  best  methods  to  our  knowledge  regarding  speed  and  accuracy,  i.e., 
Jenkins/Traub  program  and  the  eigenvalue  method. 

Furthermore  we  give  a  simple  approach  to  improve  the  accuracy  for  spectral  factorization 
in  the  case  there  are  double  roots  on  the  unit  circle.  Finally  we  briefly  consider  the  inverse 
problem  of  root  Rnding,  i.e.,  computing  the  polynomial  coefficients  from  the  roots  which  may 
lead  to  surprisingly  large  numerical  errors. 


‘This  research  was  partially  supported  by  Alexander  von  Humboldt-Stiftung  and  by  AFOSR  under  grant  F49620- 
1-0006  funded  by  DARPA 
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1  Introduction 


Finding  polynomial  roots  rapidly  and  accurately  is  an  important  problem  in  various  areas  of  signal 
processing  such  as  spectral  factorization,  phase  unwrapping,  forming  a  cascade  of  second  order 
systems  etc.  [1,  3,  4,  16,  18].  There  exists  a  large  number  of  different  methods  for  hnding  all 
polynomial  roots  either  iteratively  or  simultaneously.  Most  of  them  yield  accurate  results  only  for 
small  degrees  or  can  treat  only  special  polynomials,  e.g.,  polynomial  with  real  roots. 

One  of  the  two  best  general  purpose  root  hnder  is  the  Jenkins/Traub  method  [7].  It  works  with 
the  polynomial  itself.  The  required  memory  is  proportional  to  0(n).  The  maximum  degree  yielding 
reasonable  accuracy  is  60-80  (see  Sec.  3). 

The  second  method  is  called  eigenvalue  method  and  works  with  the  so-called  companion  matrix 
formed  with  the  polynomial  coefficients.  The  polynomial  roots  are  the  eigenvalues  of  this  companion 
matrix  which  can  be  found  with  high  accuracy  by  use  of  the  EISPACK  routines  [17].  These  are 
the  best  known  programs  for  solving  general  eigenvalue  problems.  The  required  memory  and 
computation  time  are  proportional  to  O(n^)  and  O(n^),  respectively.  As  Toh  and  Trefethen  point 
out  in  [19]  one  can  only  hope  to  get  no  worse  conditioned  problem  than  the  underlying  roothnding 
problem  by  using  the  eigenvalue  approach.  This  means  one  should  solve  the  polynomial  zerohnding 
problem  and  not  the  eigenvalue  problem  if  the  interesting  parameters  are  the  roots  of  a  given 
polynomial. 

We  present  a  method  for  hnding  all  polynomial  roots  of  an  arbitrary  complex  valued  polynomial. 
It  turns  out  to  be  faster  and  at  least  as  accurate  as  the  best  known  methods  for  nearly  all  polyno¬ 
mials  we  have  used  for  testing.  It  basically  consists  of  a  combination  of  two  well-known  iterative 
methods,  i.e.,  Muller’s  and  Newton’s  method  [14].  The  hrst  is  numerically  robust  and  yields  an 
estimate  for  the  root  working  with  the  actual  debated  polynomial.  In  a  second  step  this  root  is 
rehned  using  the  quadratically  converging  Newton’s  method  for  the  original  polynomial.  We  do  not 
assume  any  special  structure  of  the  polynomial  and  take  no  extra  precautions  for  multiple  roots. 
The  latter  should  be  done  by  the  user  who  can  exploit  additional  knowledge  about  the  roots.  As  an 
example  we  give  a  straightforward  approach  for  spectral  factorization  where  the  problem  of  double 
roots  on  the  unit  circle  may  occur 

The  paper  is  organized  as  follows.  In  Sec.  2  we  describe  Muller’s  and  Newton’s  method  and  some 
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steps  to  stabilize  the  implemented  program.  The  efficiency  and  reliability  of  all  three  methods  are 
compared  in  Sec.  3  where  we  use  several  test  polynomials  with  known  roots.  We  discuss  how  to 
deal  with  double  roots  on  the  unit  circle  in  the  case  of  spectral  factorization  in  Sec.  4.  There  we 
additionally  consider  the  problem  of  multiplying  the  roots  to  get  the  original  polynomial.  This 
simple  looking  procedure  may  lead  to  completely  perturbed  polynomial  coefficients  even  in  well 
conditioned  cases.  The  main  results  are  summarized  in  Sec.  5 

2  Description  of  the  Procedure 

We  consider  a  complex  valued  polynomial 

n  n 

P(x^  —  ^  ^  —  Pn  Pn  7”  ^  (^) 

i/=0  i/=l 

of  degree  n.  The  problem  we  address  is  to  hud  all  n  roots  Sj,  of  P(x)  as  accurately  and  fast  as 
possible.  We  work  with  a  given  and  hxed  computation  accuracy,  i.e.,  IEEE-P754-floating  point 
standard  (accuracy  ~  2.2  •  10“^®).  The  relative  accuracy  of  the  resulting  roots  has  to  be  compared 
to  this  number.  Eor  the  special  case  of  real  coefficients  pj,,  complex  roots  Sj,  must  appear  as  complex 
conjugate  pairs. 

We  have  chosen  a  combination  of  Muller’s  and  Newton’s  method  since  both  of  them  can  be  used 
to  hud  complex  roots.  In  the  following  we  give  a  pseudo  code  of  our  program: 

1.  Check  polynomial  and  return  if  erroneous  input. 

2.  If  polynomial  degree  is  1  or  2  ^  compute  roots;  return. 

3.  CallmonicO. 

4.  While  degree  of  deflated  polynomial  >  2. 

•  Call  muller  0 . 

•  Call  newtonO . 

•  Call  poldeflO . 

5.  Compute  root(s)  of  deflated  polynomial. 
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6.  Call  newt on() . 


In  a  first  step  the  coefficients  are  formally  checked,  e.g.,  possible  roots  at  zero  are  determined  and 
debated,  leading  zeros  are  cancelled  etc.  In  the  case  of  a  first  or  second  order  polynomial  the 
well-known  explicit  formulae  are  used.  Only  for  degree  ra  >  2  we  choose  our  iterative  procedure.  At 
first  the  routine  monicO  yields  a  polynomial  with  =  1.  Then  the  routine  mullerO  computes 
an  estimate  for  a  root  of  the  actual,  debated  polynomial  Pdefiix)  which  contains  all  roots  of  P{x) 
but  the  roots  found  up  to  the  actual  iteration  step  {k). 

In  a  second  step  the  estimate  resulting  from  mullerO  is  used  as  the  initial  value  of  Newton’s 
method.  It  is  simple  and  known  to  have  at  least  quadratical  convergence  near  the  solution  (for  single 
roots).  We  use  the  original  polynomial  to  avoid  errors  introduced  by  the  defiation  process.  Once 
Newton’s  method  has  converged  the  resulting  root  is  debated  (and  possibly  its  complex  conjugated 
for  real  valued  polynomials).  This  defiation  procedure  is  repeated  until  the  resulting  polynomial  is 
of  degree  two  or  one.  An  estimate  for  its  root(s)  can  be  obtained  by  using  the  well-known  explicit 
formula  and  is  refined  by  Newton’s  method. 

2.1  Muller’s  method 

We  have  chosen  Muller’s  method  for  computing  an  initial  estimate  of  the  root  because  it  has  the 
following  two  properties:  One  is  a  good  convergence  to  a  reasonable  estimate  of  a  root.  The  second 
is  the  possibility  to  get  complex  roots  even  when  initialized  with  real  values  in  opposite  to  other 
methods,  e.g.,  Newton’s  method.  The  convergence  speed  is  super  linear  (1.84  for  single  roots). 
General  convergence  for  this  method  has  not  been  proven. 

Muller’s  method  extends  the  idea  of  the  secant  method  which  works  with  a  linear  polynomial  to 
a  quadratical  polynomial.  Given  three  previous  estimates  and  x^^'>  for  an  unknown 

root  we  compute  a  new  value  by  determining  one  of  the  roots  of  a  parabola  P(x)  which  interpolates 
P(x)  in  these  three  “old”  points.  This  is  illustrated  by  Fig.  1.  The  corresponding  iteration  formulae 
are  [14] 

hk  = 

—  hfc/hfc— 1 

=  rfcP(a:W)  -  rk(l  +  rk)P(x(’^-^y)  +  rlPix^^-^)) 


Tk 

Ak 
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(2) 

(3) 

(4) 


Figure  1:  Iteration  step  of  Muller’s  method. 

Bk  =  (2rfc  +  l)P(a;W)-(l  +  rfc)2p(a;('=-i))  +  r^P(a;('=-2))  (5) 

Ck  =  (l  +  rfc)P(a;W)  (6) 

-2Ck 

%  =  - , 

Bk  ±  ^Bl  -  AAkCk 

a;P+i)  =  :,{k)^hkqk.  (8) 

The  new  estimate  is  determined  such  that  it  is  the  one  root  of  P{x)  closer  to  An 

example  is  given  in  Fig.  1  for  k  =  2.  In  the  case  the  denominator  vanishes  qk  is  chosen  as  \qk\  =  1 
with  arbitrary  phase.  The  algorithm  is  initialized  with  x^^'^  =  1,  x^^'^  =  —1,  and  =  0. 

For  the  practical  convergence  of  this  method  we  have  to  include  additional  measures  which  are 
summarized  in  the  following  pseudo  code  of  our  program  implementing  Muller’s  method: 

1.  Call  initializeO . 

2.  Repeat  twice. 

(a)  While  iteration  counter  <  ITERMAX  and  noise  counter  <  NOISEMAX  and  root  not 
found. 


•  Call  root_of _parabola() . 

•  Call  iteration_equation() . 
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•  Call  compute_function() . 

•  Call  check_x_value() . 

(b)  Call  root_check() . 

(c)  If  root  good  enough  return. 

After  initializing,  the  main  part  is  repeated  twice  with  different  starting  values  if  the  hrst  result  is 
not  good  enough.  The  main  loop  stops  whenever  one  of  the  following  criteria  is  fulhlled:  First  we 
give  a  maximum  number  of  iterations  considering  the  case  of  very  slow  convergence.  Second  we 
stop  when  the  iteration  is  dominated  by  noise.  This  means  that  we  get  only  minor  improvements 
in  the  range  of  the  computer  accuracy  during  a  hxed  number  of  successive  iterations.  Of  course 
the  program  stops  when 
^.(fc+l)  _  ^(k) 
x(k+k) 

holds  (root_check() ),  where  e  is  some  small  number  depending  on  the  computer  accuracy. 

The  hrst  step  of  the  main  loop  is  computing  the  roots  of  the  parabola  Eqs.  (2)-(7).  This  is 
followed  by  the  iteration  equation  (8).  We  have  observed  that  values  qk  computed  according  to 
Eq.  (7)  may  yield  too  large  changes  of  which  possibly  leads  to  another  root  and  causes  slow 
convergence.  This  can  be  circumvented  (and  is  actually  implemented  in  our  program)  by  allowing 
a  hxed  maximum  relative  increase  of  \qk\  from  one  iteration  step  to  the  next. 

Before  the  new  function  value  is  evaluated  we  estimate  to  avoid  overhow.  This  is  done 

by  checking  n  ■  log;^Q  If  an  estimate  indicates  a  value  greater  than  the  maximum  possible 

computer  number  we  choose 

+  (10) 

instead  of  Eq.  (8)  and  repeat  this  until  no  overhow  occurs.  This  means  we  go  back  closer  and  closer 
to  the  old  value 

In  a  last  step  the  actual  value  |P(a;(^+^))|  is  compared  to  the  best  value  until  the  current  iteration 
step  and  it  substitutes  the  latter  if  it  is  smaller. 

The  stopping  criteria  e,  ITERMAX,  . . .  were  determined  on  an  experimental  basis  to  get  a 
program  which  is  reliably  and  fast. 


<  e 


(9) 
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2.2  Newton’s  method 


Newton’s  method  is  well-known  and  works  with  the  following  simple  iteration  formula 


x(k+l)  =  2.(fc)  _  ^  ^  (k)  ^ 

’  p'ixi^yy 


(11) 


which  is  illustrated  by  Fig.  2. 


Figure  2:  Iteration  step  of  Newton’s  method. 


It  yields  the  root  of  the  tangent  through  the  point  P(a;(^)))  of  the  previous  iteration  step 

and  is  initialized  with  the  result  of  Muller’s  method.  Similar  to  this,  additional  measures  must  be 
included  to  improve  the  performance  of  the  program.  This  can  be  seen  by  the  following  pseudo 
code  of  our  implementation: 

1.  While  iteration  counter  <  ITERMAX. 

•  Call  fdvalueO . 

•  If  |P(a;(^+l))|  <  \P(x^,n)\  Xrmn  =  . 

•  If  |P'(a;(^+^))|  y  0  and  |P(a;(^+^))|/|P'(a;(^+^))|  < 

else  Aa;(^)  = 
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•  If  |Aa;(^)|/|a;(^)|  <  e  or  noise  counter  >  NOISEMAX. 

o  If  in  the  case  of  a  real  valued  polynomial  Im(a;)  <  S  ^  Im(a;mm)  =  0;  return. 


2.  If  in  the  case  of  a  real  valued  polynomial  Im(a;)  <  S  ^  Im(a;mm)  =  0. 


As  in  Muller’s  method  we  give  a  maximum  number  of  iterations.  The  hrst  step  of  the  main  loop 
is  computing  the  function  value  and  its  derivative  at  the  actual  point  This  substitutes  the 

minimum  value  up  to  the  actual  iteration  step  x^in  if  it  yields  a  smaller  function  value.  We  do  not 
permit  too  large  changes  by  using  the  new  improvement  according  to  Eq.  (11)  only  if  it  is 

smaller  than  the  old  one  of  the  previous  iteration  step.  Otherwise  the  latter  one  is  used  to  avoid 
that  the  algorithm  switches  to  another  root. 

The  algorithm  is  stopped  when  it  is  dominated  by  noise  (see  Muller’s  method)  or  when 

Aa;(fc) 

e  =  - 

^min 


where  is  that  leading  to  the  minimum  |P(a;)|  up  to  the  actual  iteration  step  and  e 

again  depends  on  the  computer  accuracy.  In  the  case  of  a  real  valued  polynomial  (i.e.,  all  pj,  are 
real)  for  every  complex  root  also  its  complex  conjugate  is  a  root  which  can  be  deflated  together. 
Consequently  we  have  to  decide  whether  a  root  with  a  very  small  imaginary  part  is  to  be  seen  as  a 
real  or  a  complex  root.  We  assume  to  have  a  real  root  if  the  imaginary  part  is  less  than  S  which  we 
have  chosen  as  half  of  the  computer  accuracy.  To  avoid  this  decision  which  may  lead  to  errors  one 
can  simply  multiply  the  real  valued  polynomial  by  the  imaginary  unit.  In  this  case  the  program 
interpretes  the  polynomial  as  a  complex  valued  polynomial  and  it  consequently  deflates  only  one 
root  at  each  iteration  step. 

It  is  interesting  that  the  quotient  e  introduced  in  (12)  can  be  used  as  an  estimate  of  the  relative 
accuracy  of  the  actual  root.  Our  program  yields  the  corresponding  maximum  value  of  all  computed 
roots  which  is  discussed  in  the  following  section. 
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3  Test  of  the  New  Algorithm 

3.1  General  Considerations 

After  the  short  description  of  our  implementation  we  have  to  prove  the  efficiency.  We  point  out 
again  that  we  do  not  try  to  construct  a  new  algorithm  with  nice  theoretical  properties  but  a 
tool  which  works  reliable  and  fast  in  many  practical  applications.  Our  approach  to  verify  the 
performance  is  oriented  at  the  ideas  of  Jenkins  and  Traub  [10].  They  propose  to  chose  a  lot  of  test 
polynomials  with  known  roots  testing  programs  for  different  weaknesses.  A  root  hnder  is  the  better 
the  smaller  the  difference  between  the  correct  roots  and  the  determined  ones.  We  use  a  normalized 
version  of  this  criterion 

Xij  ^vrain 

e  =  max  - 

where  is  the  value  computed  by  the  root  hnder  and  Sj,  is  the  corresponding  exact  value. 

This  number  is  computed  for  every  polynomial.  Additionally  we  determined  the  necessary  CPU 
time  on  an  HP  Apollo  workstation  9000/705.  We  compare  the  results  of  our  program  regarding 
to  speed  and  accuracy  with  two  of  the  best  root  hnder  programs  to  our  knowledge.  These  are  the 
Jenkins/Traub  program  [7]  and  the  eigenvalue  method  based  on  EISPACK  [17]  in  the  version  of 
MATE  AH. 

We  do  not  show  results  of  factoring  actual  transfer  functions  since  we  do  not  know  the  correct 
roots  and  cannot  give  an  objective  measure  for  the  accuracy.  However,  we  successfully  computed 
the  roots  of  many  FIR  hlters.  To  give  an  example  our  program  determined  all  roots  of  a  degree 
1000  FIR  low-pass  hlter  within  8.35s  with  an  estimated  error  e  =  1.8  •  10“^®.  The  computed  roots 
in  the  stopband  which  should  have  magnitude  one  have  a  maximum  distance  from  the  unit  circle 
of  1.11  •  10“^®.  That  means  they  are  exact  within  computer  accuracy. 

3.2  Test  polynomials 

The  following  polynomials  used  for  the  detection  of  different  properties  of  a  root  hnder  were  pro¬ 
posed  by  Jenkins  and  Traub  [10].  We  use  several  polynomials  for  each  property  but  we  leave 
out  polynomials  with  random  coefficients  for  the  same  reasons  we  did  not  examine  the  factoriza¬ 
tion  of  transfer  functions.  We  conjecture  that  in  this  case  one  gets  similar  results  to  those  of  the 
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last  polynomials  we  present  in  this  paper  (equi distantly  distributed  roots  on  the  unit  circle)  since 
polynomials  with  random  coefficients  tend  to  have  a  similar  root  distribution  [2,  18]. 

3.2.1  Check  of  the  Stopping  Criterion 

The  hrst  polynomial 

Pi{x)  =  B(x  —  A)(x  +  A)(x  —  1)  (14) 

shows  the  effect  of  very  large  or  small  roots  (A)  and  very  large  or  small  polynomial  coefficients 
(B),  respectively  on  the  stopping  criterion.  Table  1  shows  the  result  for  all  three  methods,  where 
eE,  ej,  ejv  means  the  error  of  the  eigenvalue  method,  the  Jenkins/Traub  method,  and  our  new 
method  according  to  Eq.  (13).  e  means  the  corresponding  estimate  computed  by  our  program  and 


A 

B 

^E 

ej 

ejv 

e 

tE 

tj 

tN 

1  • 10“^° 

1 • 10“^° 

0 

3.009-10-2° 

0 

0 

0.01 

0.05 

0.02 

1 • 10“^° 

1  •  10^° 

0 

3.009-10-2° 

0 

0 

0.01 

0.02 

0.01 

1  •  10^° 

1 • 10-1° 

0 

1.177-10-° 

0 

3.584-10-1^ 

0.02 

0.02 

0.01 

1  •  10^° 

1  •  10i° 

0 

-1.177-10-° 

0 

3.584-10-1^ 

0.01 

0.03 

0.01 

Table  1:  P\{x) 

=  B( 

'x  —  A){x  +  A){x  —  1). 

tE,  tj,  tN,  are  the  necessary  CPU  times  in  seconds.  As  can  be  seen  our  program  computes  all  roots 
exactly  just  as  the  eigenvalue  method  but  in  opposite  to  the  Jenkins/Traub  method  which  yields 
relatively  large  errors  for  large  values  of  A. 

Furthermore  we  use  the  polynomial 

n 

P2{x)=X{{x (15) 

i/=0 

with  more  and  more  zeros  close  to  0  for  larger  values  n.  The  results  for  ra  =  5,  7  are  summarized 
in  Table  2.  Again  our  program  yields  the  most  accurate  results  but  this  time  together  with  the 
Jenkins/Traub  method.  The  large  value  fjv  can  be  explained  by  the  fact  that  in  the  case  n  =  7  the 
main  loop  of  Muller’s  method  has  to  be  computed  twice. 

From  the  test  polynomials  above  it  can  be  concluded  that  our  method  has  the  fewest  problems 
concerning  the  stopping  criterion. 
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n 

eE 

ej 

ejv 

e 

tE 

tj 

tN 

5 

6.939-10-1® 

1.735  -  10-1® 

1.735-10-1® 

1.307-10-1® 

0.01 

0.01 

0.03 

7 

1.388-10-1® 

1.735  -  10-1® 

1.735-10-1® 

1.576  -  10-1® 

0.02 

0.02 

0.10 

n 

Table  2:  P2{x)  =  H  “  10^"). 

i/=0 


3.2.2  Check  of  Convergence 

For  the  check  of  problems  concerning  the  convergence  we  choose  the  polynomial 

n 

P?,{x)  =  W{x  -  v)  (16) 

v=l 

which  has  a  surprisingly  ill  conditioned  numerical  behaviour  (cf.  Wilkinson  [20]).  The  results  are 
depicted  in  Table  3.  We  make  the  following  observations:  The  accuracy  of  all  three  methods  is 
comparable.  For  degrees  greater  than  8  our  method  always  yields  the  best  results.  The  estimate 
e  is  close  to  the  actual  error.  For  larger  degrees  our  method  takes  the  most  CPU  time  since  the 
main  loop  in  Muller’s  method  is  computed  twice. 

3.2.3  Multiple  or  Clustered  Roots 

Multiple  roots  or  roots  close  to  each  other  lead  to  numerically  ill  conditioned  polynomials.  Nearly 
every  root  hnder  has  difficulties  to  compute  these  with  high  accuracy.  As  Wilkinson  points  out 
in  [20]  the  limiting  accuracy  of  a  root  with  multiplicity  m  is  where  e  means  the  computer 

accuracy.  This  means  in  our  case  (accuracy  ~  2.2  •  10“^®)  that  a  double  root  can  be  determined 
up  to  an  accuracy  of  about  10“®  as  long  as  no  special  methods  are  used  in  this  case  (cf.  Sec.  4.1). 
Furthermore  this  means  without  additional  knowledge  it  cannot  be  decided  whether  two  roots 
differing  by  about  10“®  are  two  separated  roots  or  incorrectly  determined  double  roots.  However, 
as  Wilkinson  points  out  in  [20]  (and  is  conhrmed  by  our  experience)  this  must  not  affect  the 
accuracy  of  well  conditioned  roots  even  for  a  polynomial  with  several  multiple  roots  (see  below). 
We  have  used  the  following  polynomials. 

P^ix)  =  (a;  -  0.1)®(a;  -  0.5)(a;  -  0.6)(a;  -  0.7) 

P^{x)  =  {x-0.1)‘^{x-0.2f{x-0.3f{x-0A) 

Pq(x)  =  (a;  -  0.1)(a;  -  1.001)(a;  -  0.998)(a;  -  0.99999) 


12 


(17) 

(18) 
(19) 


n 

eE 

ej 

ejv 

e 

tE 

tj 

tN 

1 

0.000 

0.000 

0.000 

2.220  - 

10- 

-16 

0.01 

0.01 

0.00 

2 

0.000 

2.220  • 

10- 

-16 

0.000 

2.220  - 

10- 

-16 

0.02 

0.01 

0.00 

3 

5.921  • 

10- 

-16 

6.518  • 

10- 

-20 

2.220- 

10- 

-16 

4.534- 

10- 

-25 

0.02 

0.02 

0.00 

4 

1.021  • 

10- 

-14 

1.136  • 

10- 

-19 

1.776- 

10- 

-15 

1.776  - 

10- 

-15 

0.02 

0.02 

0.01 

5 

5.003  • 

10- 

-14 

1.510  • 

10- 

-14 

1.554- 

10- 

-15 

1.776  - 

10- 

-14 

0.01 

0.02 

0.02 

6 

2.236  • 

10- 

-13 

4.269  • 

10- 

-14 

5.695- 

10- 

-14 

1.089  - 

10- 

-13 

0.02 

0.02 

0.01 

7 

8.777- 

10- 

-13 

1.305  • 

10- 

-13 

6.370- 

10- 

-13 

2.425  - 

10- 

-13 

0.01 

0.02 

0.02 

8 

1.180  • 

10- 

-11 

5.340  • 

10- 

-13 

1.217- 

10- 

-12 

1.929  - 

10- 

-12 

0.02 

0.02 

0.03 

9 

1.062  • 

10- 

-10 

9.705  • 

10- 

-12 

5.760- 

10- 

-12 

7.370  - 

10- 

-12 

0.03 

0.03 

0.04 

10 

4.365  • 

10- 

-11 

8.977- 

10- 

-11 

1.604- 

10- 

-11 

4.312  - 

10- 

-11 

0.02 

0.02 

0.05 

11 

5.448  • 

10- 

-10 

3.513  - 

10- 

-10 

1.363- 

10- 

-10 

1.772  - 

10- 

-10 

0.02 

0.02 

0.07 

12 

6.751 

•  10- 

-9 

8.332  - 

10- 

-10 

1.976- 

10- 

-10 

2.205  - 

10- 

-9 

0.02 

0.02 

0.08 

13 

1.121 

•  10- 

-7 

1.397 

-  10- 

-8 

4.252 

-  10- 

-9 

4.507- 

10- 

-9 

0.03 

0.03 

0.09 

14 

4.421 

•  10- 

-8 

4.147 

-  10- 

-8 

1.372 

-  10- 

-8 

5.249  - 

10- 

-8 

0.03 

0.02 

0.12 

15 

3.848 

•  10- 

-7 

2.725 

-  10- 

-7 

9.540 

-  10- 

-8 

5.714- 

10- 

-8 

0.03 

0.03 

0.13 

n 

Table  3:  P^ix)  =  H  “  ^)- 

i/=i 


P^{x)  =  {x-  0.001)(a;  -  0.01)(a;  -  0.1)(a;  -  0.1  +  A-i)(x  -  0.1  -  A-i)(x  -  l)(x  -  10)  (20) 
The  results  are  depicted  in  Table  4  and  5. 

Again  the  accuracy  of  all  three  methods  is  comparable  where  the  Jenkins/Traub  method  has  ad¬ 
vantages  for  Prix).  The  eigenvalue  method  always  yields  the  worst  results.  The  resulting  accuracy 
lies  in  the  expected  range  of  the  limiting  accuracy. 

We  want  to  remark  that  although  Steiglitz  and  Dickinson  conjecture  the  polynomial  P(x)  = 
(^.loo  _  _  should  stop  any  root  solver  in  its  tracks.”  our  program  is  able  to  hnd  all  roots 

of  P(x)  =  —  l)®(a;  -|-  0.5)(a;  —  0.5)(a;  —  2)(x  —  3)  with  a  maximum  error  of  4  •  10“^  (for  the 

multiple  roots).  All  single  roots  could  be  computed  up  to  computer  accuracy  which  conhrms  the 
proposition  above  about  the  effect  of  multiple  roots. 
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u 

eE 

ej 

ejv 

e 

tE 

tj 

tN 

4 

2.124-  10“® 

3.435-10“® 

8.912  -  10“® 

3.248  -  10“® 

0.01 

0.04 

0.04 

5 

7.076-  10“'^ 

1.441  -  10“® 

5.626  -  10“'^ 

1.492  -  10““^ 

0.03 

0.03 

0.09 

6 

1.529-10“® 

2.918-10“'^ 

1.002  -  10“® 

4.116  -  10“® 

0.01 

0.03 

0.04 

Table  4:  Polynomials  Pj^{x),  v  =  4,5,6. 


A 

eE 

ej 

ejv 

e 

tE 

tj 

tN 

0 

1.627- 

10- 

-5 

4.430  - 

10“® 

6.978  - 

10“® 

2.406- 

10“® 

0.01 

0.03 

0.04 

1  - 

10- 

-10 

1.627- 

10- 

-5 

4.384- 

10“® 

6.978  - 

10“® 

2.406- 

10“® 

0.02 

0.01 

0.04 

1 

-  10- 

-9 

1.364- 

10- 

-5 

4.063  - 

10“® 

8.620  - 

10“® 

1.841  - 

10“® 

0.02 

0.02 

0.04 

1 

-  10- 

-8 

1.610  - 

10- 

-5 

1.013  - 

10“" 

4.988  - 

10“® 

1.145- 

10“® 

0.02 

0.02 

0.04 

1 

-  10- 

-7 

9.612  - 

10- 

-6 

9.902  - 

10“" 

4.863  - 

10“® 

3.632- 

10“® 

0.02 

0.04 

0.03 

1 

-  10- 

-6 

2.407- 

10- 

-5 

9.989  - 

10“® 

2.918  - 

10“® 

1.289- 

10“® 

0.03 

0.02 

0.04 

Table  5:  Polynomial  Prix). 


3.2.4  Stability  of  Deflation 

Stability  of  the  deflation  process  means  that  the  roots  of  the  deflated  polynomial  are  close  to  those 
of  the  original.  Since  we  rehne  all  roots  using  the  original  polynomial  our  method  is  expected  to 
yield  good  results.  We  choose  the  polynomials 


Psix) 

=  (x  —  A)ix  —  l){x 

-I' 

(21) 

M-l 

3M 

Pgix) 

=  n  ) 

■  n  “  0.9  -  e^^). 

(22) 

v=M 

where  the  latter  has  a  distribution  of  roots  similar  to  the  transfer  function  of  an  FIR  lowpass 
hlter.  Table  6  and  7  show  the  results  and  Fig.  3  depicts  the  accuracy  and  the  CPU  time  for  P^ix) 
depending  on  the  degree  n. 
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A 

eE 

ej 

ejv 

e 

tE 

tj 

tN 

1  •  10^ 

2.274-  10-1® 

2.483-10-1® 

0.000 

1.139-10-1® 

0.01 

0.02 

0.00 

1  •  10® 

2.328  •  10-1® 

2.220-10-1® 

2.118  -  10-1® 

1.164-10-1® 

0.02 

0.01 

0.00 

1  •  10® 

4.441  •  10-1® 

2.068-10-1® 

2.068  -  10-1® 

1.192-10-1® 

0.01 

0.02 

0.00 

Table  6:  Psix)  =  (x  —  A)(x  —  l)(x  —  ^). 


M 

eE 

ej 

ejv 

e 

tE 

tj 

tN 

5 

2.357-10-1® 

1.959  -  10-1® 

2.758  -  10-1® 

1.820-10-1® 

0.16 

0.03 

0.05 

10 

3.903-10-1® 

2.227-10-1® 

3.084-10-1® 

1.825-10-1® 

0.92 

0.05 

0.15 

12 

4.242-10-1® 

1.771  -  10-1 

4.934-10-1® 

1.364-10-1® 

1.43 

0.07 

0.19 

15 

6.690-10-1® 

1.552-10-1 

8.723  -  10-1® 

9.812-10-1^ 

2.63 

0.11 

0.27 

20 

2.534-10-11 

2.396-10-1 

1.061  -  10-1® 

1.825-10-1® 

6.27 

0.18 

0.39 

25 

1.300-10-1® 

1.613-10-1 

4.939  -  10-1® 

1.875-10-1® 

12.18 

0.30 

0.51 

50 

1.803  -  10-® 

3.510-10® 

2.481  -  10-1® 

1.121  -  10-11 

95.34 

26.11 

2.31 

Table  7:  Pgix 

II 

1  1— 1 1 

1 

e'2M j .  n 

v=M 

-  0.9-ej 

V7T 

2M  ) 

Figure  3:  Results  of  Jenkins/Traub  (  +  ),  eigenvalue  (o),  and  our  method  (x)  for  Pdix).  (a)  Actual 
and  estimated  accuracy  e  and  e  (—  •  — );  (b)  CPU  time. 
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Our  method  always  yields  the  most  accurate  results.  It  is  better  than  a  factor  of  1000  for 
Pg(x),  n  =  200  (M  =  200)  and  furthermore  much  faster  for  increasing  degree  up  to  a  factor  of 
45  compared  to  the  eigenvalue  method.  It  is  remarkable  that  the  Jenkins/Traub  program  yields 
much  less  accurate  results  even  for  Pg(x),  ra  =  40  (M  =  10)  and  completely  useless  values  for  larger 
degrees.  As  we  will  see  this  is  typical  for  the  program. 

3.2.5  High  Order,  Well  Conditioned  Polynomials 

As  a  last  example  we  consider  high  degree  polynomials  up  to  ra  =  10000  for  the  very  well  conditioned 


polynomials 

Piq(x)  =  a;”  —  1 

(23) 

Pii(x)  =  a;”  -|-  1. 

(24) 

The  results  are  summarized  in  Table  8  and  9  and  graphically  depicted  in  Figs.  4  and  5. 


n 

eE 

ej 

ejv 

e 

tE 

tj 

tN 

10 

8.083  •  10-1® 

2.428-10-1® 

2.289  -  10-1® 

1.579  -  10-1® 

0.04 

0.02 

0.02 

20 

1.226  •  10-1® 

5.431  -  10-1® 

5.979  -  10-1® 

1.247-10-1® 

0.08 

0.06 

0.04 

50 

2.109  •  10-1® 

6.785-10-12 

1.144-10-1® 

2.220  -  10-1® 

0.34 

0.11 

0.11 

70 

2.047-10-1® 

3.049  -  10-1 

6.280  -  10-1® 

5.155  -  10-1^ 

0.72 

0.24 

0.17 

100 

2.559  •  10-1® 

9.769  -  10-1 

1.024-10-1® 

1.868  -  10-1® 

1.77 

0.43 

0.26 

200 

3.443  •  10-1® 

- 

1.180  -  10-1® 

2.220  -  10-1® 

11.88 

- 

0.67 

500 

3.360  •  10-1® 

- 

8.968  -  10-1® 

2.069  -  10-1® 

315.22 

- 

2.59 

1000 

- 

- 

1.024-10-1® 

2.202  -  10-1® 

- 

- 

8.81 

2000 

- 

- 

1.106  -  10-1® 

2.059  -  10-1® 

- 

- 

30.96 

10000 

- 

- 

1.047-10-1® 

2.196  -  10-1® 

- 

- 

944.09 

Table  8:  Pio(x)  =  x^  —  1 

Again  the  accuracy  of  the  Jenkins/Traub  program  drastically  decreases  for  relatively  small  de¬ 
grees  (n  =  50).  It  cannot  be  used  for  degrees  ra  >  60  . .  .70.  On  the  other  hand  our  method  yields 
the  best  results  with  nearly  computer  accuracy  up  to  the  degree  n  =  10000.  The  eigenvalue  method 
is  slightly  worse  regarding  the  accuracy  but  could  be  used  only  for  degrees  up  to  about  500.  This  is 
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n 

eE 

ej 

ejv 

e 

tE 

tj 

tN 

10 

7.022  •  10-1® 

1.466-10-1® 

5.661  -  10-1® 

1.784-10-1® 

0.04 

0.03 

0.02 

20 

1.422  •  10-1® 

7.462-10-1® 

1.159  -  10-1® 

1.610  -  10-1® 

0.05 

0.05 

0.04 

50 

1.355  •  10-1® 

1.082  -  10-® 

1.024-10-1® 

1.964-10-1® 

0.34 

0.14 

0.11 

70 

1.490  •  10-1® 

3.569  -  10-1 

6.280  -  10-1® 

1.110  -  10-1® 

0.76 

0.25 

0.18 

100 

2.253  •  10-1® 

3.680  -  10-1 

1.180  -  10-1® 

7.390  -  10-1^ 

1.74 

0.45 

0.32 

200 

2.811  •  10-1® 

- 

1.180  -  10-1® 

2.073  -  10-1® 

11.77 

- 

0.67 

500 

3.581  •  10-1® 

- 

1.024-10-1® 

2.059  -  10-1® 

316.12 

- 

2.64 

1000 

- 

- 

1.106  -  10-1® 

2.144-10-1® 

- 

- 

9.04 

2000 

- 

- 

1.043  -  10-1® 

1.794-10-1® 

- 

- 

31.48 

10000 

- 

- 

1.024-10-1® 

2.219  -  10-1® 

- 

- 

1036.80 

Table  9:  P\\{x)  =  a;”  +  1 


because  of  the  necessary  large  memory  to  store  the  companion  matrix  rP  ■  8Byte  (double  precision) 
which  is  2  •  10®Byte  for  n  =  500  and  8  •  10®Byte  for  n  =  10000.  The  CPU  time  for  our  method  is 
always  the  smallest  as  can  be  seen  in  Figs.  4(b)  and  5(b).  Especially  for  n  =  500  it  is  2.59s  {Piq{x)) 
and  2.64s  (Pii(a;))  compared  to  315.2s  {Piq{x))  and  316.1s  (Pii(a;)).  As  in  all  examples  above  the 
estimate  e  is  close  to  the  actual  accuracy  of  our  method. 
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Figure  4:  Results  of  Jenkins/Traub  (  +  ),  eigenvalue  (o),  and  our  method  (x)  for  Pio(x).  (a)  Actual 
and  estimated  accuracy  e  and  e  (—  •  — );  (b)  CPU  time. 


P{x)  =  x*n  + 1 


Figure  5:  Results  of  Jenkins/Traub  (  +  ),  eigenvalue  (o),  and  our  method  (x)  for  Pii(x).  (a)  Actual 
and  estimated  accuracy  e  and  e  (—  •  — );  (b)  CPU  time. 
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4  Discussion 


4.1  Spectral  Factorization 

The  problem  of  spectral  factorization  is  to  find  all  roots  of  a  symmetric  polynomial  H{x)  (a  poly¬ 
nomial  where  the  existence  of  a  root  at  Sj,  implies  the  existence  of  a  root  at  (a;*)“^)  with  the 
additional  property  that  all  roots  on  the  unit  circle  have  even  multiplicity.  In  the  following  we 
assume  they  have  multiplicity  two.  After  hnding  the  roots  of  H{x)  one  is  interested  to  form  a 
minimum  phase  polynomial  or  in  general  a  polynomial  P{x)  of  degree  n  such  that 

H{x)  =  P{x)-x'^P*{{x*)-^)  (25) 

holds. 

As  we  mentioned  earlier,  in  general  the  limiting  accuracy  for  double  roots  is  half  the  computer 
accuracy.  However,  in  the  special  case  where  we  know  their  modulus  it  is  possible  to  use  a  simple 
procedure  to  compute  them  up  to  computer  accuracy.  This  can  be  done  by  the  following  procedure 
where  we  assume  that  all  roots  on  the  unit  circle  are  well  separated  which  holds  in  nearly  all 
practical  cases. 

In  a  hrst  step  we  compute  all  zeros  of  the  original  polynomial  p[{x)  and  separate  them  into  those 
lying  inside,  outside,  and  on  the  unit  circle.  To  separate  these  three  regions  we  chose  an  annulus 
with  thickness  100  times  the  expected  accuracy  where  we  use  the  estimate  e  of  our  program.  In 
a  second  step  we  compute  the  roots  of  H'(x).  It  has  the  same  roots  on  the  unit  circle  as  the 
original  polynomial  H(x)  but  with  multiplicity  1.  Consequently,  they  can  be  computed  with  higher 
accuracy.  For  this  polynomial  we  are  only  interested  in  the  roots  lying  on  the  unit  circle  where  we 
now  chose  a  much  thinner  separating  annulus  corresponding  to  the  higher  accuracy.  If  the  number 
of  roots  of  H(x)  on  the  unit  circle  is  twice  the  number  of  those  of  H'(x)  these  can  be  immediately 
used  as  an  improved  estimate. 

We  have  implemented  this  approach  into  a  MATLAB  hie  (see  Appendix)  and  give  an  example 
by  using  the  polynomial  P(x)  =  P^ix)  with  degree  n  =  100.  We  computed  H(x)  according  to  Eq. 
(25)  and  the  corresponding  minimum  phase  part  according  to  the  method  above.  The  maximum 
error  of  the  roots  of  the  resulting  polynomial  compared  to  the  original  one  is  1.4  •  10“^“^,  which  is 
hardly  worse  than  the  achieved  accuracy  of  the  original  polynomial  (cf.  Table  7,  M  =  25).  The 
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roots  of  H{x)  on  the  unit  circle  could  be  determined  only  with  accuracy  6  •  10“®  so  that  there  is 
an  improvement  of  a  factor  5  •  10®. 

We  point  out  that  the  operation  of  computing  the  polynomial  coefficients  from  the  roots,  we  call 
it  coefficient  hnding,  must  be  done  with  care.  Otherwise  the  errors  introduced  by  it  can  be  much 
larger  than  the  given  accuracy.  This  is  investigated  in  the  following  section. 

4.2  Coefficient  Finding 

As  we  have  seen  coefficient  hnding  is  a  necessary  step  when  we  do  spectral  factorization.  The  fact 
that  computing  the  polynomial  coefficients  from  the  roots  can  lead  to  very  large  errors  although 
the  polynomial  has  the  best  numerical  condition  one  can  think  of,  seems  not  to  be  well-known. 
The  only  related  reference  we  could  hnd  is  [13].  We  want  to  give  an  intuition  of  what  problems 
may  arise  and  how  these  can  be  overcome. 

We  again  consider  the  polynomial  Pio(x)  with  all  its  roots  //  =  1 . .  .ra  on  the  unit  circle. 

Let  us  assume  we  compute  the  polynomial  coefficients  in  a  straightforward  manner  using  the  order 
implied  by  the  roots  given  above.  The  resulting  coefficients  corresponding  to  the  terms  x  . .  .x^~^ 
are  expected  to  be  zero  within  computer  accuracy.  However,  for  n  =  20,50, 100,200  the  maximum 
values  of  these  coefficients  are  8.7  •  10“®^,  1.6  •  10“®,  8  •  10^,  and  5.5  •  10®°®.  This  at  the  hrst  glance 
surprising  result  can  be  explained  as  follows. 

In  the  process  of  computing  the  polynomial  coefficients  we  have  intermediate  polynomials  with 
roots  only  on  a  sector  of  the  unit  circle.  These  are  known  to  be  ill  conditioned  [11].  Furthermore 
they  have  coefficients  with  a  large  dynamical  range.  This  last  property  leads  to  errors  in  the 
next  step  where  the  intermediate  polynomial  has  to  be  convolved  with  a  hrst  order  polynomial  and 
numbers  of  different  orders  have  to  be  added.  The  last  observation  shows  an  easy  way  to  circumvent 
this  problem.  One  merely  has  to  sort  the  roots  in  such  a  way  that  each  intermediate  polynomial 
has  approximately  equidistantly  distributed  roots.  These  lead  to  coefficients  with  a  small  dynamic 
range  of  the  coefficients. 

If  ra  is  a  power  of  two  this  can  be  easily  done  by  interpreting  the  indices  //  =  1 ...  ra  of  the  roots  as 
binary  numbers,  reverse  the  order  of  the  digits  and  interpret  this  again  as  a  decimal  number.  This 
procedure  is  known  as  bit  reversal  or  van  der  Corput  sequence  [15,  6].  In  the  case  n  is  not  a  power 
of  two  we  simply  continue  the  sequence  1 . .  .ra  to  the  next  larger  power  of  two.  Then  we  proceed  as 
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described  and  cancel  all  values  larger  than  n  from  the  resulting  sequence.  As  an  example  consider 
the  vector  [1 . .  .20]  which  becomes 

[1  17  9  5  13  3  19  11  7  15  2  18  10  6  14  4  20  12  8  16]. 

Using  this  “reindexing”,  the  error  for  n  =  20,50,100,200  is  10“^®,  4  •  10“^®,  6  •  10“^®,  3  •  10“^“^, 
which  is  a  drastical  improvement  compared  to  the  results  above. 

Although  this  bit  reversal  is  primarily  suited  for  equidistantly  distributed  roots  on  the  unit  circle 
it  can  also  be  successfully  used  for  coefficient  hnding  of  an  FIR  hlter  since  these  often  have  a 
similar  root  distribution.  However,  there  is  a  generalization  called  Leja  ordering  for  arbitrary  root 
locations  [13]  which  is  computationally  more  expensive  but  yields  slightly  improved  results  for  FIR 
hlters  compared  to  the  bit  reversal.  The  corresponding  MATLAB  functions  can  be  found  in  the 
Appendix. 

5  Conclusion 

The  following  conclusions  can  be  drawn  from  the  examples  given  in  Sec.  3:  All  three  methods  yield 
comparable  results  regarding  speed  and  accuracy  when  working  with  “difficult”  low  degree  (n  <  20) 
polynomials.  For  these  the  CPU  time  is  on  a  low  and  comparable  level.  The  accuracy  of  our  method 
is  always  better  than  that  of  the  eigenvalue  method  and  better  than  that  of  the  Jenkins/Traub 
method  in  most  cases.  For  larger  degrees  (ra  >  30  . .  .40)  the  accuracy  of  the  Jenkins/Traub  method 
drastically  decreases  and  it  yields  useless  results  for  n  >  60... 70.  The  accuracy  of  our  method 
is  better  than  that  of  the  eigenvalue  method  in  every  case  sometimes  by  more  than  a  factor  of 
1000  (cf.  Pg(x)).  Furthermore  the  CPU  time  of  our  method  increases  much  slower  than  that  of  the 
eigenvalue  method,  e.g.,  it  is  faster  than  a  factor  of  hundred  for  n  =  500.  This  fact  and  the  linear 
increasing  memory  (depending  on  n)  compared  to  a  quadratical  needed  for  the  eigenvalue  method 
makes  it  possible  to  hud  roots  in  considerable  time  even  for  high  degree  polynomials  (~  1000s  for 
a  degree  10000  polynomial). 

To  summarize  the  results,  our  program  seems  to  have  no  drawbacks  (comparable  to  good  methods 
for  low  degrees)  but  essential  advantages  regarding  accuracy  and  speed  which  is  especially  true  for 
for  large  degrees.  A  C  version  of  the  program  can  be  obtained  by  the  authors. 
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Furthermore  we  gave  a  powerful  approach  for  spectral  factorization  which  considerably  improves 
the  accuracy.  Finally  we  have  considered  the  inverse  problem  to  factorization,  i.e.,  hnding  the 
polynomial  coefficients  from  the  roots.  We  showed  that  large  errors  can  result  from  this  process 
and  we  gave  a  simple  method  to  minimize  them  nearly  up  to  computer  accuracy. 


A  M-File  for  Spectral  Factorization 


function  [p_out]  =  factorize(p_in, string) 

%  function  [p_out]  =  factorize (p_in , string) 

% 

%  Input:  p_in,  string 

I 

%  Output :  p_out 

% 


% 

% 

% 

% 

% 

% 

% 

% 

% 

% 


Description:  Program  yields  (spectral)  factorization  of  the  polynomial 
p_in  regarding  to  the  the  unit  circle.  For  reliable 
results  p_in  should  have  single  roots  off  the  unit  circle 
and  double  roots  on  the  unit  circle.  If  the  optional  input 
variable  string  contains  at  least  one  'a'  the  maximum  phase 
portion  is  determined  (taking  the  double  roots  of  p_in  on  the 
unit  circle  once  +  all  roots  of  p_in  outside  the  unit  circle)  . 
In  all  other  cases  the  function  yields  the  minimum  phase 
portion . 


%  subroutines:  ransort.m  bitrev.m  isodd.m  rootsl.mex 


'/oFile  Name:  factorize.m 

'/oLast  Modification  Date:  XG'/o  XU'/o 

'/oCurrent  Version:  "LWL 

'/File  Creation  Date:  Tue  Sep  7  09:29:06  1993 
/Author:  Markus  Lang  <lang@dsp . rice . edu> 

% 

/Copyright:  All  software,  documentation,  and  related  files  in  this  distribution 
/  are  Copyright  (c)  1993  Rice  University 

/ 

/Permission  is  granted  for  use  and  non-profit  distribution  providing  that  this 
/notice  be  clearly  maintained.  The  right  to  distribute  any  portion  for  profit 
/or  as  part  of  any  commercial  product  is  specifically  reserved  for  the  author. 

/ 

/Change  History: 

/ 


p_in  =  p_in( : ) ' ; 
n  =  length(p_in)  -  1; 
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%  find  roots  of  original  and  differentiated  polynomial 
[r_orig,e]  =  rootsl (p_in) ;  r_orig  =  r_orig(:).';  error_orig  =  e; 

[r_dif,  error_dif]  =  rootsl ( (n : -1 : 1) . *p_in(n+l : -1 : 2) ) ;  r_dif  =  r_dif 

%  find  roots  on,  inside  and  outside  the  unit  circle  (original  pol.) 

%  and  sort  by  angle 

ind_orig_u  =  f ind(abs (1-abs (r_orig) ) <100*error_orig) ; 
ind_orig_i  =  f ind(abs (r_orig) <abs (l-100*error_orig) ) ; 
ind_orig_o  =  f ind(abs (r_orig) >abs (l+100*error_orig) ) ; 
r_orig_u  =  r_orig(ind_orig_u) ;  [dum,ind]  =  sort (angle (r_orig_u) ) ; 
r_orig_u  =  r_orig_u(ind) ;  n_orig_u  =  length(r_orig_u) ; 
r_orig_i  =  r_orig(ind_orig_i) ;  [dum,ind]  =  sort (angle (r_orig_i) ) ; 
r_orig_i  =  r_orig_i (ind) ;  n_orig_i  =  length(r_orig_i) ; 
r_orig_o  =  r_orig(ind_orig_o) ;  [dum,ind]  =  sort (angle (r_orig_o) ) ; 
r_orig_o  =  r_orig_o (ind) ;  n_orig_o  =  length(r_orig_o) ; 

%  check  multiplicity  of  roots  on  unit  circle 
if  isodd(n_orig_u) ; 

disp('There  are  roots  on  unit  circle  with  odd  multiplicity! !'); 
return 

end 

%  check  whether  all  roots  could  be  sorted 
if  n_orig_u+n_orig_i+n_orig_o  ~=  n; 

disp('Not  all  roots  could  be  sorted'); 
return 

end 

%  find  roots  on  the  unit  circle  (differentiated  pol.)  and  sort  by  angle 
ind_dif_u  =  f ind(abs (1-abs (r_dif ) ) <100*error_dif ) ; 
r_dif_u  =  r_dif (ind_dif _u) ;  [dum,ind]  =  sort (angle (r_dif_u) ) ; 
r_dif_u  =  r_dif_u(ind) ;  n_dif_u  =  length(r_dif _u) ; 

%  check  whether  the  numbers  of  roots  on  the  unit  circle  of  both  polynomials 
%  fit  together 
if  2*n_dif_u  ~=  n_orig_u 

disp('The  numbers  of  roots  on  unit  circle  do  not  fit  together! !'); 
return 

end 

%  check  whether  maximum  or  minimum  phase  polynomial  is  desired 
r_out  =  [r_dif_u  r_orig_i] ; 
if  exist (' string ' )  ==  1; 

if  length(f ind(string==' a' ) )  >  0; 

r_out  =  [r_dif _u  r_orig_o] ; 
end 

end 
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%  compute  output  polynomial 

[dum,ind]  =  sort (angle (r_out) ) ;  r_out  =  r_out(ind); 
p_out  =  poly (r_out (ransort (length(r_out) ) ) ) ; 


function  odd  =  isodd(x) 

%  function  odd  =  isodd(x) 

I 

%  function  yields  a  matrix  odd  which  is  1  for  every  odd 
%  value  of  abs(x)  and  0  else. 

I 

%  ml,  20.8.1992 

% 

%  Copyright  Lehrstuhl  fuer  Nachrichtentechnik  Erlangen,  ERG 
%  e-mail:  int@nt.e-technik.uni-erlangen.de 

[m,n]  =  size(x) ; 
odd  =  zeros(m,n); 

ind  =  f ind( (x-1) /2==f ix( (x-1) /2) ) ; 
odd(ind)  =  1  +  odd(ind) ; 

B  M-File  for  Bit  Reversal 

function  [ind]  =  ransort (n) 

%  function  [ind]  =  ransort (n) 

% 

%  function  computes  an  index  vector  ind  of  length  n  which  contains 
%  all  integers  l:n  but  in  a  quasi  random  order. 

% 

%  m-file  bitrev  must  be  available 
%  author:  m.  lang  17.12.91 

%  This  is  very  advantageous  for  example  for  the  function  poly: 

%  compute  z  =  roots ([1  zeros (l:m)  1]);  p  =  poly(z) . 

%  There  is  a  large  error  in  p,  which  stems  from  the  second  operation 
%  poly.  If  the  roots  z  are  sorted  with  increasing  angle  in  zs,  then 
%  ph  =  poly (zs (ransort (n) ) )  yields  much  less  error  than  p. 

I 

%  ransort  uses  bitreversal  which  is  known  as  van  der  corput  sequence 
%  in  the  theory  of  uniform  distributed  sequences. 

m  =  ceil (log2 (n) ) ; 
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ind  =  1 : 2"m; 

ind  =  bitrev(ind); 

ind  =  ind(f ind(ind<=n) ) ; 


function  re  =  bitrev(rf eld) 

% 

%  BITREV  Y  =  BITREV(X)  returns  the  vector  X  in  bitreversed  order 

I 

I 

%  Author  :  Raimund  Meyer  15.12.87 

%  Revised:  Rai  26.08.88 

%  Copyright  (C)  1987-1991  Lehrstuhl  fuer  Nachrichtentechnik, 

%  University  of  Erlangen,  ERG 

dim  =  length(rf eld) ; 
diml  =  dim  -  1 ; 

j  =  1; 

for  i  =  l:diml 
if  i  <  j 

xt  =  rf eld( j ) ; 

rfeld(j)  =  rfeld(i) ; 
rfeld(i)  =  xt ; 
else 
end 

k  =  dim/2; 
while  k  <  j 
j  =  j  -  k; 
k  =  k/2; 

end 

j  =  j  +  k; 

end 

re  =  rfeld; 

C  M-File  for  Leja  Ordering 

function  [x_out]  =  leja(x_in) 

%  function  [x_out]  =  leja(x_in) 

I 

%  Input:  x_in 

% 

%  Output :  x_out 

I 

%  Program  orders  the  values  x_in  (supposed  to  be  the  roots  of  a 

%  polynomial)  in  such  a  way  that  computing  the  polynomial  coefficients 
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%  by  using  the  m-file  poly  yields  most  accurate  results. 

%  Try,  e.g., 

%  z=exp(j*(l:100)*2*pi/100) ; 

%  and  compute 

%  pi  =  poly(z); 

%  p2  =  poly(leja(z)) ; 

%  which  both  should  lead  to  the  polynomial  x"100-l.  You  will  be 

%  surprised! 

I 

%  ref.:  Nachtigal,  Reichel,  Trefethen.  "A  Hybrid  GMRES  Algorithm  for 

%  Nonsymmetric  Linear  Systems.  SIAM  J.  Matr.  Anal,  and  Appl .  ,  13, 

%  796-825,  July  1992. 

'/oFile  Marne:  leja.m 

'/oLast  Modification  Date:  XG'/o  XU'/o 

'/oCurrent  Version:  '/oM'/o 

'/oFile  Creation  Date:  Mon  Mov  8  09:53:56  1993 
'/oAuthor:  Markus  Lang  <lang@dsp .  rice  .  edu> 

I 

'/oCopyright :  All  software,  documentation,  and  related  files  in  this  distribution 
%  are  Copyright  (c)  1993  Rice  University 

I 

'/oPermission  is  granted  for  use  and  non-profit  distribution  providing  that  this 
'/onotice  be  clearly  maintained.  The  right  to  distribute  any  portion  for  profit 
'/oOr  as  part  of  any  commercial  product  is  specifically  reserved  for  the  author. 

% 

'/oChange  History: 

I 


X  =  x_in(:).';  n  =  length(x); 

a  =  x(ones (1 ,n+l) , : ) ; 
a(l, :)  =  abs(a(l, :))  ; 

[duml,ind]  =  max(a(l , 1 :n) )  ; 

dum2  =  a(:,l);  a(:,l)  =  a(:,ind);  a(:,ind)  =  dum2 ; 

x_out(l)  =  a(n,ind); 

a(2,2:n)  =  abs (a(2 , 2 : n) -x_out (1) ) ; 

for  l=2:n-l 

[duml,ind]  =  max(prod(a(l : 1 , 1 : n) ) ) ;  ind  =  ind+1-1; 
if  l~=ind 

dum2  =  a(:,l);  a(:,l)  =  a(:,ind);  a(:,ind)  =  dum2 ; 

end 

x_out(l)  =  a(n,l); 

a(l+l , (1+1) : n)  =  abs (a(l+l , (1+1) : n) -x_out (1) ) ; 
end 

x_out  =  a(n+l ,  :  )  ; 


26 


References 


[1]  Aliphas,  Amnon,  S.  Shankar  Narayan,  and  Allen  M.  Peterson.  Finding  the  zeros  of 
linear  phase  hr  frequency  sampling  hlters.  IEEE  Transactions  on  Acoustics,  Speech, 
and  Signal  Processing,  31:729-734,  June  1983. 

[2]  Arnold,  G.  fiber  die  Nullstellenverteilung  zufalliger  Polynome.  Mathematische 
Zeitschrift,  92:12-18,  1966. 

[3]  Chen  Xiangkun,  and  Thomas  W.  Parks.  Design  of  optimal  minimum  phase  hr  hlters 
by  direct  factorization.  EURASIP  Signal  Processing,  10:369-383,  1986. 

[4]  Daubechies,  Ingrid.  Ten  Tectures  on  Wavelets.  SIAM,  Philadelphia,  PA,  1992. 

[5]  Frenzel,  Bernhard.  Unter.suchung  von  Algorithmen  zur  Bestimmung  von  Polynomnull- 
stellen  (Examination  of  Algorithms  for  Einding  Polynomial  Roots,  in  German).  Studi- 
enarbeit  (senior  project).  Institute  of  Communication  Theory,  University  of  Frlangen, 
Frlangen,  1993. 

[6]  Hlawka,  F.  The  Theory  of  Uniform  Distribution.  AB  Academic  Publisher,  Berkhamsted, 
1984. 

[7]  Jenkins,  M.  A.  Algorithm  493  zeros  of  a  real  polynomial.  ACM  Transactions  on  Math¬ 
ematical  Software,  1:178-,  June  1975. 

[8]  Jenkins,  M.  A.  and  J.  F.  Traub.  A  three-stage  algorithm  for  real  polynomials  using 
quadratic  iteration.  SIAM  Journal  on  Numerical  Analysis,  7:545-566,  1970. 

[9]  Jenkins,  M.  A.  and  J.  F.  Traub.  Zeros  of  a  complex  polynomial.  Communications  of 
the  ACM,  15:97-99,  February  1972. 

[10]  Jenkins,  M.  A.  and  J.  F.  Traub.  Principles  of  testing  polynomial  zerohnding  programs. 
ACM  Transactions  on  Mathematical  Software,  1:26-34,  March  1975. 


27 


[11]  Lang,  Markus.  Ein  Beitrag  zur  Phasenapproximation  mit  Allpdssen  (Phase  Approxi¬ 
mation  hy  Allpass  Filters,  in  German).  PhD  thesis,  University  of  Erlangen-Niirnberg, 
Erlangen,  Germany,  1993. 

[12]  Muller,  D.  E.  A  method  for  solving  algebraic  equations  using  an  automatic  computer. 
Math.  Tables  Aids  Comput,  10:208-215,  1956. 

[13]  Nachtigal,  Noel  M.,  Lothar  Reichel,  and  Lloyd  N.  Trefethen.  A  hybrid  GMRES  algo¬ 
rithm  for  nonsymmetric  linear  systems.  SIAM  Journal  on  Matrix  Analysis  and  Appli¬ 
cations,  13:796-825,  July  1992. 

[14]  Press,  William  H.  et  ah  Numerical  Recipes  in  C.  Cambridge  University  Press,  Cam¬ 
bridge,  1992. 

[15]  Reichel,  Lothar,  and  Gerhard  Opfer.  Chebyshev- Vandermonde  systems.  AMS  Mathe¬ 
matics  of  Computation,  57:703-721,  October  1991. 

[16]  Schmidt,  C.  E.  and  L.  R.  Rabiner.  A  study  of  techniques  for  Ending  the  zeros  of  linear 
phase  hr  digital  filters.  IEEE  Transactions  on  Acoustics,  Speech,  and  Signal  Processing, 
25:96-98,  Eebruary  1977. 

[17]  Smith,  B.  T.  et  ah  Matrix  Eigensystem  Routines  —  EISPACK  Guide,  volume  6  of 
Tecture  Notes  in  Computer  Science.  Springer,  1976. 

[18]  Steiglitz,  Kenneth,  and  Bradley  Dickinson.  Phase  unwrapping  by  factorization.  IEEE 
Transactions  on  Acoustics,  Speech,  and  Signal  Processing,  30:984-991,  December  1982. 

[19]  Toh,  Kim-Chuan,  and  Lloyd  N.  Trefethen.  Pseudozeros  of  polynomials  and  pseudospec¬ 
tra  of  companion  matrices.  Technical  Report  93-1360,  Department  of  Computer  Science, 
Cornell  University,  Ithaca,  N.  Y.,  1993. 

[20]  Wilkinson,  J.  H.  Rounding  Errors  in  Algebraic  Processes.  Prentice-Hall,  Englewood 
Cliffs,  N.  J.,  1963. 


28 


